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ABSTRACT 

Formulas  for  calculating  Bessel  functions  of  integral  order  and 
complex  argument  are  derived  in  this  report.  Calculations  based  on 
these  formulas  are  not  subject  to  the  loss  of  significant  figures 
which  occurs  in  the  Taylor  and  Neumann  series  when  the  argument  is 
large  and  the  order  is  small. 

To  calculate  J  (z),  select  an  integer  m  >  n  and  m  I:1  such 

that  |.J  (zlj  <<  1,1  (z)|.  Calculate  J  (z)  and  ,J  (•■)  from  the  lavlor 
m  1  n  1  m  m+ 1 

series,  then  calculate  J  (z)  from  the  recurrence  relation.  A  simi’ar 

n 

procedure  is  used  to  calculate  I^fz). 

To  calculate  Kn(z),  express  the  quotient  Qn(z)  -  jlzj/K^fz 

in  terms  of  two  Gauss  continued  fractions.  The  individual  functions 

K  (z)  and  K  Az,  are  obtained  from  Q  (z)  and  the  Wronskian  relation 
n  n  - 1  xn 

i  volving  K  ( z  j ,  K  ,(z),  1  t  z)  ,  and  I  ,/.zj.  A  similar  nroeedure 
*  nv  ’  n-1  n  n-1 

involving  llankel  functions  is  used  to  calculate  Y^fz)  and  j ( z ) . 
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Cancellation  of  significant  figures  is  very  severe  in  the  Neumann 
series  for  when  n  is  small  and  2  is  large  and  real.  At  the 

present  time,  this  difficulty  is  generally  overcome  by  using  multiple 
precision  arithmetic,  quadrature  formulas,  or  continued  fractions  with 
very  involved  terms.  The  Gauss  continued  fractions  used  in  this  report 
are  simple  in  form,  rapidly  convergent,  and  are  not  subject  to  excess¬ 
ive  round-off  error. 


TABLE  OF  CONTENTS 


Page 

ABSTRACT . .  3 

LIST  OF  SYMBOLS . . . 

I.  INTRODUCTION .  (J 

II.  ORDINARY  BF.SSEL  FUNCTIONS  OF  THE  FIRST  KIND .  10 

III.  MODIFIED  BESSEL  FUNCTIONS  OF  THE  FIRST  KIND .  12 

IV.  MODIFIED  BESSEL  FUNCTIONS  OF  THE  SECOND  KIND . IS 

V.  ORDINARY  BESSEL  FUNCTIONS  OF  THE  SECOND  KIND . ID 

VI.  RESULTS  AND  CONCLUSIONS . 21 

REFERENCES .  21 

APPENDIX . 25 

DISTRIBUTION  LIST . 5' 


S 


!  1ST  OF  SYMBOLS* 


>’ 


k ,f  ,m,n 


Jn(z) 

YnU) 

V21 

K  fz) 
n 

(2) 

H  l  J  (:.l 
n 

f (a ,b ,uj 

,mUJ 

Qn(Z) 

Pnj(.-i 

ln(2)ui 


I:j  ( a ,  b ; ii ) 
Oj (a,b;u) 


0 


rfz) 


t 


,  •  ,  i  0 

complex  variable,  z  =  x  +  ,y  -  pc 

■“cal  part  of  z 
imaginary  part  of  z 
absolute  value  of  z 
argument  of  z 
integers 

ordinary  Bessel  function  of  the  first  kind 

ordinary  Bessel  function  of  the  second  kind 

modified  Bessel  function  of  the  first  kind 

modified  Bessel  function  of  the  second  kind 

llankel  functions,  sometimes  called 
Bessel  functions  of  the  third  kind 

KullS  form  of  the  confluent  hypergeometric  function 

Whittaker's  form  i  the  confluent  hypergeometric  fun  tier, 

a  quotient.  Q  (z)  =  K  ,  (z)/K  (z) 
n  n  n- 1  n 

quotients  of  Hankel  functions; 
see  equations  (4b),  (47) 

a  complex  variable  defined  in  terms  of  z;  see 
equation  (*?1) 

(iauss  continued  fractions.  See  equotio.s  (30)-i  , 

of  the  order  of 
the  gamma  function 
Filler's  constant 


*:he  rotation  of  Chapter  ‘J ,  Ref.  5,  is  used  whenever  practical. 
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a 

normalizing  factor  in  Miller's  algorithm 

V2) 

Ln(z)  =  aJn(z)  +  6Vn(z),  Eq.  (5) 

e 

a  small  quantity 

Re 

the  real  part  of 

Im 

the  imaginary  part  of 

R{F(z) } 

the  real  part  of  F(z),  Appendix  A 

I{F(z)} 

the  imaginary  part  of  F(z),  Appendix  A 

F(z) 

a  function  of  z 

w 

dependent  variable  in  the  differential  equation 
,2 

d  w  „ ,  ,  A 

— =-  -  F(z)w  =  0 

dz 

I 
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I .  INTRODUCTION 


Linear  boundary  value  problems  in  cylindrical  coordinates  may 

frequently  be  solved  in  terms  of  Bessel  functions  of  integral  order. 

Certain  problems  in  elasticity,  viscoelasticity,  and  fluid  flow  require 

Bessel  functions  of  a  complex  argument  in  order  to  satisfy  the  boundary 

conditions  in  a  rigorous  manner.  Taylor  and  Neumann  series  are 

generally  used  to  calculate  these  functions  when  the  argument  is  small 

1  2  3* 

or  moderate  in  size.  ’  ’  However,  severe  cancellation  occurs  unless 
all  the  terms  of  the  series  have  the  same  sign.  Consequently,  multiple 
precision  arithmetic  is  required  in  order  to  obtain  final  results  of 
the  required  accuracy. 

Cancellation  is  a  severe  form  of  round-off  error  which  occurs  if 
the  sum  of  a  series  is  small  compared  with  the  largest  term.  Several 
significant  figures  may  be  lost  by  subtraction  if  the  calculations  are 
carried  out  with  a  fixed  number  of  significant  figures.  This  diffi¬ 
culty  is  illustrated  by  the  alternating  series 


-x  ,  121  3 

e  =  1  -  X  ♦  y  X  -  jr  x  +  •  •  • 

since  e  X  is  small  when  x  is  large.  It  is  evident  that  cancellation 
does  not  occur  when  evaluating  the  series 


1  . 

X  +  T  X 


1 


Hence,  to  obtain  e  X,  we  first  calculate  o  and  then  obtain  e  X  by 
d i v i s i on . 

Miller  has  shown  that  cancellation  does  not  occur  if  stable 
recurrence  relations  are  used  to  calculate  Bessel  functions  of  the 
first  Kind.  -\n  inulogous  procedure  cannot  be  used  for  calculating 
Bessel  functions  of  the  second  Kind,  .»s  the  corresponding 
relations  are  unstable.”* 


R  c  i  c  ic  h  c  1 2  at  c  t  is  i  e  d  c »:  pa  a  c  2  4 . 

Afvjot-ifhr  di-ScuiSed  cn  pajoi  555  and  59',  Kc5.  5. 
Stability  d  vxiicui  vcutvncc  tct'afto:  o  a\  jcujsocf  c:  pa,?c  Kill  c; 
ti:c  h\t-.cdu:  ti.cn,  Re  i .  5. 


In  this  paper,  the  quotient  of  two  modified  Ek  _,el  functions  of 
the  second  kind  is  expressed  in  terms  of  two  Gauss  continued  fractions. ^ 
The  original  function  is  then  calculated  from  the  Wronskian  relation 
and  the  functions  of  the  first  kind.  This  procedure  eliminates  cancel¬ 
lation,  so  tl.at  accurate  values  of  the  function  can  be  obtained  for  a 
wide  range  of  order  and  argument.  A  similar  procedure  is  used  to  calcu¬ 
late  ordinary  Bessel  functions  of  the  second  kind. 

If  we  require  a  sequence  of  functions  J^(z),  J^+^(z)  ...  Jn+^(z)> 

we  calculate  J  ■  (z)  and  J  ,(z),  then  calculate  the  functions  of 
n  +  t  n  +  i+  l 

lower  order  from  the  recurrence  relation.  Similarly,  if  Y  (z), 

1  n 

Y  , ( z )  ...  Y  . (z)  are  required,  we  calculate  Y  (z)  and  Y  ,(z),  then 
n-t-1  n+i  1  n  n+1 

calculate  functions  of  higher  order  from  the  recurrence  relations.  The 
routine  being  developed  at  the  Ballistic  Research  Laboratories  (BRL1 
will  provide  the  pair  of  Bessel  functions  required  to  start  the 
recurrence  process. 


II.  ORDINARY  BESSEL  FUNCTIONS  OF  THE  FIRST  KIND 
These  functions  can  be  calculated  accurately  from  the  series* 


J 

n 


It) 


n 


1 ! (n+1) 


:: (n+ i j (n+:) 


3  !  (n  +  1)  (n  +2'1  (n  +  3) 


U) 


if  :i;  is  small,  or  if  n  is  very  large  compared  with  | z ;  .  However,  if 
z  is  real  and  large  compared  with  n,  the  n*^  term  of  the  series  is 
large  compared  with  J  (z) ;  consequent »y ,  severe  cancellation  occurs  in 
summing  the  series. 


Ke  will  describe  Miller's  Algorithm  briefly,  as  it  is  the  basis 
of  most  subsequent  work  in  the  area.  Ke  assume 

m  =  o 


m+ 


=  i 

m 


U’) 

13) 


iovcc.  cCr.J  -ccu.\  tcx afconA  „tv  c 

5  <1*1  ci  7. 


:  otvio'.t'u  stem  tc ao*ovcc j 
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where  m  •  >  n,  m  >  >  !z|.  The  functions  L  ,  (z),  L  .  (zi  ...  L  (  .)  are 

m- 1  m-d  n 

generated  from  the  recursion  formula 


‘t'-k-iW  *  =  2kLk(=> 


(4) 


We  express  L  (z)  as  a  linear  combination  of  J  ( z J  and  Y  (z): 
1  nr  mv  nv  ’ 


I.  (z)  =  aJ  (z)  +  bY  (z) 
m  m  m 


(5) 


Then,  since  J  (z) .  Y  (z)  ,  and  I.  (z)  all  satisfy  the  same  linear  recur- 
m  m  m  1 

rence  relation,  we  have 


LnU)  =  aJnU)  +  6YnUj 


(6) 


When  |  z  j  '■  k,  J  ^  ( - )  increases  and  Y^(z)  decreases  as  k  decreases;  con¬ 
sequently  the  second  term  in  Hq.  (6)  is  negligible,  and  we  have 


L n  1  z  j  'v  a.)  ( :  J 


(') 


The  normalizing  factor  i  is  obtained  from  the  Neumann  series  of  an 
elementarv  function. 


k-n 


or 


fin 


k  =  n 


k  =  n 


k  k 


1S1 


It  is  assumed  that  I  a  l  l;.)  is  negligible.  !'h x s  procedure  is  very 

et loot ive  provided  cancel lat ion  does  no’  occur  in  summing  the  series  in 
kq.  (81.  However,  it  is  difficult  to  ‘itui  a  single  function  f(n,z) 
which  sitisties  this  reou  >  remen  t  when  noth  n  ■•ad  :  u  v  widely. 


Ip.  contrast  with.  Miller’s  algorithm,  the  author  uses  iq.il'  to 
calculate  d  (z  l  and  J  ,  ( z  1 ,  where  m  •  •  z  and  m  n .  The  terms  of 
the  series  are  moderate  in  size,  ana  cancellation  is  negligible.  The 
function  ’  (:)  is  calculated  bv  repeated  use  of  the  recurrence  relation 

n  * 


z[Jk.1(z)  ♦  J};  +  1(z^]  =  2kJk(z), 


(9) 


.-.irting  with  k  -  n. 

The  integer  n  was  obtained  from  the  polynomial  approximation 


.1  ,2  ,1  ,4  ,1  ,21 

^  (■)  p)  (9  p)  (•>  b  i  / 1  a-j 

1  *  P  (1  _ 2  _  l _  _ ^ _ 1  (10) 

l,r.  “  2^,  ’  li(n+l)  2 !  (n+1)  (n+T)"  ’’’  £!’(n+l) .  .  .  (n +1)  ’ 


where  p  =  |zj  and  i  is  the  smallest  integer  for  which 


fl  ,1 
(t  0} 

4 

rr(n^~7(^rf 


(id 


Then 


m  =  n  +  t 


(12) 


The  small  constant  t  was  adjusted  by  trial.  The  final  value  was  chosen 

just  small  enough  to  guarantee-  the  accuracy  of  J  (z)  in  the  range  of 

n 

interest  * 


III.  HOD  I  FI. fcD  SESSEL  FUNCTIONS  OF  THE  FIRST  KIND 
A  similar  procedure  is  used  to  calculate  I  (z)  from  the  series 


In(z)  = 


,n  1 
2  n! 


{1 


1  2 

(T  z  ) 


rl 

cj  i.  ) 


1 !  (n+1)  +  2! (n+1)  (n+2) 


,1  ,b 

( -  7  1 

_ ]2  22 _ _ 

3:  (n+l)’(n+2) (n+2) 


}  (A3) 


Obviously.,  no  cancellation  occurs  if  z  is  reel  and  positive,  as  the 
terms  of  the  series  are  all  positive.  If,  however,  z  is  a  puio 
inu.5inary  number,  severe  cancellation  will  occur  when  jzj  is  large  com¬ 
pared  with  n. 

We  calculate  I ^ ( z )  and  I  ^(z)  where  the  integer  m  is  given  by 
Eq.  (12).  The  function  *  (z)  is  calculated  by  successive  applications 
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* 


of  the  recurrence  relation 

-  -kV^  (14) 

starting  with  k  =  m.  This  recurrence  relation  is  stable  for  decreasing 
index . 

A  number  of  functions  which  have  Taylor  series  expansions  and  reci 
ronce  relations  which  are  stable  for  decreasing  index  can  be  calculated 
by  Miller's  algorithm  or  procedures  similar  to  the  one  outlined  above. 
If,  however,  the  function  has  a  logarithmic  singularity  at  the  origin, 
no  simple  method  of  calculation  based  on  series  and  recurrence  relations 
appears  to  be  available. 


IV.  MODIFIED  BESSEL  FUNCTIONS  OF  THE  SECOND  KIND 
These  functions  can  be  calculated  from  the  Neumann  series 


KqU) 


[y  -  loge  2  +  logp  =]I0(z 


(j  z)2 

___ 


.1  , 4  1  ,6 

*■2  1  2  ^  11 
-j — t (1  +  ■=■)  +  ■; — ^ — 5 — : — - — —  (1  +  +  — )  +  .  .  . 

1  * 2 *  1  * 2  v  2  l*2*3*l*2*o  2  o 


(IS) 


Kn(z)  = 


[y 


log,. 


+  log  zj 


n-1 

i  x 

r=0 


(-1)' (n-r-1) ! 


(|) 


n-2r 


n-2r 


(~ z) 


r !  (n+r) 


[1 


1  +4+...+  -J— ]  (16) 

r  2  n+rJ  v  J 


r=0 


if  j  z |  is  small,  or  if  |z|  >  n.  However,  severe  cancellation  occurs  if 
z  is  large,  real,  and  greater  than  n.  In  this  range,  we  have  the 
approximat i ons 


Kn(z) 


!nU)  - 


:wz 
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so  that  ( z)  is  small  and  I  (z)  is  large.  It  is  evident  from  Eq.  (16) 

that  K^(z)  will  be  calculated  from  the  difference  of  two  large  and 
nearly  equal  numbers. 

Although  the  Neumann  series  is  useful  when  n  is  large  compared  with 
j  z i ,  functions  of  lower  order  cannot  be  calculated  accurately  from  the 
recurrence  relation 

z[K  , (2)  -  K  ,  (z) ]  =  2nK  (zl  (17) 

n+1  n-1  J  n 

as  the  differences  of  nearly  like  numbers  also  occur  in  the  course  of 
these  calcuic ' ions .  It  appears  that  K^(z)  should  be  calculated  from  a 
formula  which  does  not  separate  the  analytic  and  logarithmic  parts  of 
the  function. 

This  requirement  is  satisfied  by  the  integral  representation* 


Kn«  *  (2f,2“TT  /  e'“  t’  *  27>”T  du  ‘18> 

I  (n*j)  J 

o 

which  is  valid  provided  z  does  not  lie  on  the  negative  half  of  the  real 
axis.  It  is  related  to  a  form  of  the  confluent  hypergeometric  function 
discussed  by  Wall.** 


f (a,b;v) 


,  / •  -u  a-1, 

1  /  e  u  du 

r<a>  J  Cl  *  vu)¥ 
0 


(19) 


We  see  that 


1 

Kn(z)  =  (27) 2  e"Z  f(r  b;v) 


(20) 


*Vage  HTf'Wfn 
**Page  £C£,  Rejj.  8 
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1 

a  =  2  +  n 


Similarly 


(21) 

(22) 

(23) 


Kn-l(z)  =  (2l)2  e'Z  -  1,  b  ♦  1;  v) 
We  define  the  quotient  function  by  the  equation 


Qn(z)  =  Kn-l(z)/Kn(z) 


(24) 


(25) 


Then 

Qn(Z)  =  f(a  -  1,  b  +  1;  v)/-'a,b;v)  (26) 

We  reduce  the  e  .  press  ion  on  the  right  to  a  form  that  can  be 
expressed  in  terms  of  Gauss  continued  fractions.  The  function  f(a,b;v) 
satisfies  the  recurrence  relations 

f(a,  b;  v)  =  f(a  +  1,  b;  v)  +  bvf(a  +  1,  b  +  1;  v)  (27) 

f(a,  b;  v)  =  f(a,  b  +  1;  v)  +  avf(a  +  1,  b  +  1;  v)  (28) 

Now  replace  a  by  a  -  1  and  b  by  b  +  1  in  Eq.  (27).  We  find 

f(a  -  1,  b  +  1;  v)  =  f(a,  b  +  1 ;  v)  +  (b  *  1)  vf(a,  b  +  2;  v)  (29) 

so  that 


n  -  f(a,  b  *■  1;  v)  (b  +  1)  vf(a,  b  2;  v) 
^nl  j  f(a,b;v)  f(a,b;v) 


On  re-arranging  this  expression,  we  find 

n  =  f(a,  b  ♦  1;  y)  (b  +  1)  vf(a,  b  *  2;  v) 

V  }  "  f(a,b;v)  f (a,  b  +  1;  v) 


(30) 


IS 


or 


where 


Qn(z)  =  Pj (a,b ;v)  1  +  v(b  +  1)  G1(a5b;v)j 


ft  r 


_  ,  v  ..  f(a,  b  +  1;  v) 

Fi(a’b;v)  = 


«  i 1 


G.(a,b;v)  *  f 
1 v  ’  '  ’  f(a,  b  +  1;  v) 

The  functions  Fj(a,b;v)  and  G^(a,b;v)  can  be  expressed  in  terras  of  Gauss 
continued  fractions.  We  have* 


f(a,b;v) 
f(a,  b  -  l';'  v) 


,  ,  (a  ♦  l)v 

,  ,  ft  » 

,  ,  (a  *  2)v 

mum 

l  +  • 


F^a.bjv) 


(b  +  l)v 


(a  +  l)v 

! ,  p+iiv 
1  +  • 


!,  RejS-  6. 


t 


anti 


Gj (a,b;v) 


(b  •»  2)v _ 

,  (a  +  i)v 

+  7  (b  +  TjT 

I  i - 

1  +  • 


To  compute  these  continued  fractions,  we  assume  that 

r  f  ,  ,_f(a+f-l,b+£-l;v) 

F£(a,b;v)  -  -f(a  +  £  _  lf  b  +  l  _  2;  v) 


G^(a,b;v) 


f (a  +  £  -  1 ,  b  +  l;  v) 
f  (a +  l  - 1 ,  b  +  t  -  \ ;  v) 


Then  we  can  prove  that 


1  +  (b  +  i)v  p£+1 (a,b;v) 

F^(a,b , v)  1  +  (b  +  Zjv  F^+1(a,b;v)  +  (a  *  i  -  l)v 


1  +  (b  +  l  +  1 ) v  G^+1  (a ,b ; v) 
G£(a,b;V)  =  I  +  (b  *  l  *  1 ) v  G^+1(a,b;v)  +  (a  *  L  -  l)v 

In  order  to  start  the  iterative  procedures,  we  assume  that 
Ff  +  i  (a>b  ;v)  =  0 
G£+1 (a,b;v)  =  0 


l  =  t 


The  integer  t  is  chosen  sufficiently  large  to  keep  the  truncation 
ms  x 

error  within  the  required  bounds. 
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rhe  individual  functions  K  (z)  and  K  , fz)  are  obtained  from  the 

n  1  n-lv 

Wronskian  relation 

♦  Viwyw  ■  {  <38> 

By  definition,  we  have 

Kn-1^  =  Kn<z>Vz> 

On  eliminating  K  . (z)  from  these  equations,  we  find 
n  ~  i 

KnU)  =  i]l  '  (z)  I  I  (z)0  (z)]'  (39) 

Convergence  of  these  continued  fractions  is  rapid  if  z  lies  in  the 
right  half  plane  and  |z|  is  large,  btr  becomes  much  slower  as  z 
approaches  the  origin  or  a  point  on  the  negative  half  of  the  real  axis. 
Moreover,  the  complex  zeros  of  K^Cz),  Kn  j(z),  and  f(a,  b  +  1;  v)  all 
lie  in  the  left  half  plane,  and  may  possibly  lead  to  division  by  zero 
for  certain  values  of  z.  Division  by  zero  cannot  occur  if  |z|  >  2,  lies 
on  the  imaginary  axis  or  in  the  right  half  plane,  and  the  integer 
^max  *s  sufficiently  large.* 

Analytic  continuation  is  v.sed  in  conjun  tion  with  the  Gauss 
continued  fraction  if 

Re  z  <  0 
and 

M  >  2 

Let  t  =  -z;  then  t  lies  in  the  righ"  half  plane,  so  that  the  functions 

K  (t)  and  K  . (t)  can  be  calculated  as  indicated  above.  K  (z)  and 
n  n-1  n 

j(z)  are  then  calculated  from  the  following  formulas  for  analytic 
continuation. 


*A  detailed  pxoo 16  given  in  Appendix  A. 
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K.(z) 

=  K  (t)  - 

ini j (t)  , 

Im  z  > 

.  o, 

j  even 

(40) 

Kj  (z) 

=  K.(t)  + 

i 71 1 j  (t) , 

Im  z  < 

o, 

j  even 

(41) 

K  (z) 

II 

1 

' _ 

»■ 

ini  .  (t.3  , 

Im  z 

I  0. 

j  odd 

(42) 

Kj(z) 

=  -Kj (t)  - 

ini  .  (t)  , 

Im  z 

<  0, 

j  even 

(43) 

These  formulas  are  obtained  by  comparing  the  Neumann  series  for  K^.(t) 
and  K.(z).  We  note  that  K^(z)  is  discontinuous  when  we  cross  the 
negative  real  axis. 


V.  ORDINARY  BESSEL  FUNCTIONS  OF  THE  SECOND  KIND 

These  functions  are  calculated  in  terms  of  the  Hankel  functions 

H  ^(z)  and  H  ^(z).  which  are  linear  combinations  of  the  ordinary 
n  n 


Bessel  functions. 

Hn(1)(z)  =  Jnfz)  ♦  iYn(z)  (44) 

hn(2)(z)  =  Jn(2)  '  iY„(z)  (45) 

The  quotient  functions 

PnU)(z)  =  llnl‘J'(z)/Hn(i)(z)  (46) 

Pn(2) (2)  =  Hnl(2)(z)/Hn(:)(z)  (47) 


are  derived  from  the  integral  representations 


I!  ll)(zl  =  (-=•) 


1  .  .  1  I  , 

£ _ ^ _ 

!’(n  ♦  4) 


/  e'“u 


1 

n— r 


1U, 


(1  *  ~)  du 


(48) 


avg  z 
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e'Uu  2(1  -  ~±)  du  (49) 

3  1 

-  —  <  aVg  Z  <  yTT 

On  referring  to  Eqs.  (19),  (21),  (22),  and  (23),  we  find 

v  =  i/2z,  Pn(1)(z)  =  i  Qn(v)  (50) 

if  Ini  ?  >  0 

v  =  -  i/2z ,  Pn(2j(z)  =  -  iQn(v)  (51) 

if  Im  z  <  0 

The  restrictions  on  Im  z  insure  that 
Re  v  _>  0 

in  both  cases.  As  Eqs.  (50)  and  (51)  taken  together  account  for  the 
entire  complex  z  plane,  additional  formulas  for  analytic  continuation 
are  not  required. 

We  now  derive  Wronskian  formulas  which  involve  only  one  type  of 
Hankel  function.  We  have 


Viw 

Vs 

)  - 

JnU) 

Yn-1< 

z)  =  -2n/z 

(52) 

On  referring 

to  Eqs. 

(44) 

and 

(45)  , 

we  find 

Vz)  « 

ilJ 

v  n 

U) 

-  (1  (1)  (z)  ] 
n 

(53) 

Vz)  . 

i[J 
k  n 

U) 

-  H  l 

n 

2)U)J 

(54) 

We  use  these 

relations  tc 

e 1 iminate 

V2> 

and  Y  ,  (z)  from  Eq.  (52) 
n  - 1 

Vi^ 

-  J 

n 

™  Hn 

(54) 

(2)  _ 


& 


-i ( z  - 


1 


r(n  .  i) 


I 
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(z)  =  -  2 n i / z 


(56) 


Jn_,(z)  Mn<^(z)  -  Jn(z)  Hn_/2» 


On  referring  to  tqs.  (50)  and  (51),  we  find 


Ha(1,(z,  . 


zl.^.^z)  -  Pnll)(z)  Jn(z)| 


H  f ! '  (z;  =  P  ( z )  li  (1^(z) 
n  n  n 


(2),=)  -  - 


:tti 


Z[J„(=)  -  Pn(2,(z)  Jn(Z)] 


(2)  ( 2 )  ( 2 ) 

Hn-1  ^  »  V  ^  V  (2J 


(57) 


(58) 

(59) 


(60) 


We  calculate  Y  ( t- )  and  ^  ( t )  from  Lqs.  (50),  (53),  (57),  and  (59), 
if  Im  :  2  0.  If  Im  z  >_0,  we  use  Lqs .  (51),  (54),  (58),  and  (60). 

This  involved  sequence  of  calculations  has  been  checked  by  extensive 
calculations  on  the  BRLHSC. 


\I.  RESULTS  AND  CONCLUSIONS 

We  have  derived  an  accurate  and  efficient  method  for  calculating 
Bessel  functions  of  the  second  kind  for  integral  order  and  complex 
argument.  I'he  formulas  given  here  have  been  programmed  in  both  the 
10RAST  and  POUT RAN  IV  programming  languages. 

In  general,  .1  (x)  and  Y  (xl  can  be  calculated  to  14  decimal  places 
■  n  n 

if  x  is  positive  and  less  than  n.  The  calculations  are  accurate  to  14 
significant  figures  if  n  •  x  .  The  functions  In(x)  and  K  (x)  can  also 
be  calculated  to  14  significant  figures.  It  is  difficult  to  check  the 
accuracy  of  the  calculations  when  the  argument  is  complex,  as  suffi¬ 
ciently  accurate  tables  are  not  available.  ITiese  remarks  apply  to  the 
range. 

2  *  x  25 
0  •  n  •  25 
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Existing  tables  were  used  whenever  possible.  In  addition,  quadrature 
formulas  based  on  Gaussian  quadrature  along  paths  of  steepest  descent 
in  the  complex  plane  were  also  derived.  The  calculations  were  time- 
consuming  but  highly  accurate.  Calculations  based  on  these  quadrature 
formulas  were  used  as  a  basis  of  comparison  whenever  tables  were  not 
avai labl e . 

In  the  routines  under  development  at  BRL,  Taylor  and  Neumann  series 
are  used  when 

No  cancellation  occurs  for  these  small  values  of  the  argument.  The 
Hankel  asymptotic  expansions  wil'  be  used  when  \z\  is  large  compared 
with  n.  The  routines  are  being  written  in  both  the  FORA? r  and  FORTRAN 
IV  programming  languages. 

The  procedure  outlined  in  fhi^  report  can  be  extended  to  Bessel 

functions  of  fractional  order,  and  also  to  Whittaker's  function  K.  (:). 

k ,  m 

The  Gauss  continued  fraction  is  most  effective  precisely  where  the 
series  development  is  plagued  with  cancellation. 
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Af'PtNDi  X  A 

A*  must  show  that  division  by  zero  will  not  occur  ii  calculating 
the  Gauss  continued  fractions  and  the  quotient  function  Q.n ( z)  provided 

j  2  j  >  2 

and 

Re  r  0. 

’the  functions  K  (z)  ana  K  ,  (:)  have  no  zeros  in  this  region.*  It 
n  n-i  6 

follows  from  Eqs .  (24)  and  (25)  that  f(a,  b;  v)  and  f (a  -  1,  b  +  1;  v) 
are  also  free  of  zeros.  It  is  shown  below  that  f(a,  b  +  1;  v)  is  also 
free  of  zeros  in  this  region;  consequently,  division  by  zeio  cannot 
occur  in  Eqs .  (26)  and  (3u>  . 

The  proof  is  based  on  an  oscillation  theorem  due  to  Einar  Hille:** 

’Met  F(z)  he  real  and  positive  when  z  is  real  avid  greater  than  x  , 
analytic  throughout  a  region  D  including  the  real  axis  for  R(z)  >  x., 
and  such  that  either 

RiF(zj)  ••  C  or  l(F(z)}  f  0 

in  0;  let  tY(z)  be  a  solution  of 

~  -  F  (z)w  -  0 

j  -  *- 
ClZ 

and  sjch  that  W(z)  ->  0  as  z  -►  in  f?  along  a  parallel  to  the  real  axis, 
then  under  very  general  assumptions,  W(z)  has  no  zero  nor  extremum 
In  D.-’ 

The  function  f(a,  b+;  v)  can  be  expressed  in  terms  of  Whittaker ’s 

function  W.  (?.},  which  satisfies  a  differential  equation  of  the 
k ,  m 

required  form.*’’* 

*Page  57 1,  fief-  ~1 . 

**Exampfe  4,  page.  527,  Ref  9;  paqt6  360  and  361,  Ref  10 . 

***ChupteA  XVI,  page;  537 -354,  Ref,  II, 


0 


(A- 1 ) 


d  *"W  (z) 

k,m 

dz*- 


l 

1  k  4 
4  +  z  + 


m 


W,  (z) 
k,mv  3 


W.  (z)  =  — 7 
k,mv  rr\ 


2  k 
e  z 


r  (t  -  k  +  m) 


1  ,1 

-  K  -  —  +  m  k  -  -  +  m 

2  ,,  u,  2 

u  (1  +  -)  e 


Now 


f(a,  b  +  1;  v)  = 


r(n  ♦  d.j 


1  3 

n  -  --  n  -  T 

2,.  .  2  -u, 

i  (i  +  vu)  e  du 


A  comparison  of  Eqs .  (A-2)  and  (A-3)  shows  that 


1  1 
V  =  7*  m  =  n  -  2*  K  =  '  2 


Eq.  (61)  now  becomes 

1 

»  H  “  — 

2*  2 _ 

dz2 


2 

d  W 


1  1  n  -  n 

4  +  2z  +  2 


(z) 


2’  n  "  2 


The  function  F(z)  of  Hille's  theorem  is 


,  1  1  n  -  n 

Ft:)  "  4  *  27  +  2 


Two  cases  arise.  If  n  =  0  or  n  =  1, 

Re{F(z)  ]  =  I  +  — 2~~ — j 
x  +  y“ 

The  function  Re{F(z)}  is  positive  outside  the  circle 
(x  +  2)2  +  yc  ~  4 


Udu  (A-2) 


(A-3) 


(A-4) 


0  (A- 5 ) 


(A-6) 
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and  negative  within  it.  The  region  D  consists  of  the  entire  complex 
z  plane  outside  the  semi-infinite  strip  hounded  by  the  semicircle 


and  the 

Wh 

F(z)  in 


and 


We  find 


of  i  f 


x  f  we  sc 

we  find 


y  =  ±  /~4  -  (x  +  2)2  ,  -2  <_  x  <  0 

lines 

v  =  ±  2,  x  <  -2. 

n  i>  >  2,  we  must  consider  both  the  real  and  imaginary  parts  of 
delineating  the  region  D.  We  assume 

Re{ F (z) }  >  0  in 

Im{F(z) }  >  0  in  D2 


lm{F(z)}  <  0  in  D7 


Re  { F  ( z ) }  =  r  + 


_X _ +  (n2  -  n)(x2  -  y2) 

4  '  2  2,  .2  2,2 

-(x  +  y  )  (x  +  y  ) 


Im{ F(z) 


} 


2(V  -  n)xy 


,  2  2,  ,  2  2,2 

2(x  +  y  )  fx  +  y  ) 


IA-7) 


(A-8) 


Ref F(z) }  >  0  if  x2  >  y‘ 


-  n) 


L 

ImF(z)  =  0 


y  =  0 


2  7 


or 

[x  +  2(n2  -  n).r  +  y2  =  4(n2  -  nj 

We  see  that  Im{F(z)}  /  0  provided  z  lies  outside  this  circle  and  does 
not  lie  on  the  x  axis.  The  regions  D^,  D^,  and  are  shown  below. 


Figure  A-l.  Zero-Free  Regions  of  W 


1 


(  =  ) 


The  region  D  contains  these  three  regions.  The  entire  right-half  plane 
and  the  imaginary  axis,  exclusive  of  the  origin,  lies  in  D. 


If  arg  z  <  r,  Whittaker's  function  has  the  asymptotic  expansion 


2  1  . 
m  -  (k  -  -) 


W.  (z)  -v  e 
k  ,mv 


1  + 


o(-4) 
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z 


(A- 9) 


and  hence  tends  to  zero  as  z  +  ®  along  a  path  parallel  to  the  x  axis. 
Consequently,  W  ^  ^(z)  has  no  zeros  in  D,  and  division  by  zero 

2’  n  ”  ~2 

will  not  occur  in  Eq.  (30).  Th '  '  conclusion  has  been  verified  by  exten¬ 
sive  calculations  on  the  BRLESC.  The  functions  occurring  in  Eq.  (33) 
and  Eq.  (34)  may  be  analyzed  in  the  same  way.  The  parameter  k  of 
Whittaker's  function  must  be  negative  or  zero  if  we  are  to  be  certain 
that  the  right  half  plane  is  free  of  zeros.  This  condition  is  satis¬ 
fied  for  the  cases  under  consideration. 

These  results  show  that  the  continued  fractions  for  F^(a,  b;  v) 
and  G^(a,  b;  v) ,  Eqs.  (31)  and  (32),  have  no  poles  in  the  right  half 
plane  or  the  imaginary  axis.  Hence,  these  continued  fractions  converge 
uniformly  in  this  region.* 
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